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ABSTRACT 

We present high-resolution computer simulations of dust dynamics and pianetesimai formation in turbulence generated by the magne- 
torotational instability. We show that the turbulent viscosity associated with magnetorotational turbulence in a non-stratified shearing 
box increases when going from 256^ to 512^ grid points in the presence of a weak imposed magnetic field, yielding a turbulent viscos- 
ity of Of ^ 0.003 at high resolution. Particles representing approximately meter-sized boulders concentrate in large-scale high-pressure 
regions in the simulation box. The appearance of zonal flows and particle concentration in pressure bumps is relatively similar at mod- 
erate (256^) and high (512^) resolution. In the moderate-resolution simulation we activate particle self-gravity at a time when there 
is little particle concentration, in contrast with previous simulations where particle self-gravity was activated during a concentration 
event. We observe that bound clumps form over the next ten orbits, with initial birth masses of a few times the dwarf planet Ceres. At 
high resolution we activate self-gravity during a particle concentration event, leading to a burst of pianetesimai formation, with clump 
masses ranging from a significant fraction of to several times the mass of Ceres. We present a new domain decomposition algorithm 
for particle-mesh schemes. Particles are spread evenly among the processors and the local gas velocity field and assigned drag forces 
are exchanged between a domain-decomposed mesh and discrete blocks of particles. We obtain good load balancing on up to 4096 
cores even in simulations where particles sediment to the mid-plane and concentrate in pressure bumps. 
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1. Introduction 

The formation of km-scale planetesimals from dust particles 
involves a complex interplay of physi cal processes, including 
most importantly collisi onal sticking (IWeid enschill ingl 11984 
; ll997l;|Dullemond & pominill2 005l). the self-gravity of the par- 
. tide mid-plane laver (ISafronov 1969; Goldreich & Ward, 1973; 
SekivaL[T998l:lYoudin & sCTl2002:.Schrapler & Henning..2004: 



Johansen et al. L l2007h . and the motion and structure of the tur- 
bulent protop lanetary disc gas ( Weiden schilling & Cuzzil 1 19931: 
iJohansen et al . 2006 : Cuzzi et al , 2008^ . 

In the initial growth stages micrometer- sized silicate 
monomers readily stick to form larger dust aggregates 
(iPoppe et all I200QI: iBlum & Wurml [2008). Further growth to- 
wards macroscopic sizes is hanipered by collisional fragmen- 
tation and bouncing (IZsom et all. 120101) . limiting the maximum 
particle size to a few cm or less (depending on th e assumed ve- 
locity t hreshold for co llisional fragmentation, see iBrauer et"all 
I2008al: iBirnstiel et al.L l2009i) . High-speed collisions between 
small i mpactors and a lar ge target constitutes a path to net 
growth (IWurm et all l2005l) . but the transport of small particles 
away from the mid-plane by tu rbulent diflTusion limits the result- 
ing growth rate dramatically (* Johansen et al , 2008). Material 
properties are also important. .Wada et al. ( 200 9) demonstrated 
efficient sticking between ice aggregates consisting of 0.1 yum 
monomers at speeds up to 50 m/s. 
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Turbulence can play a positive role for growth by co ncentrat- 
ing m m- sized particles in convection cells ([K lahr & Henningl 
1 19971) and between small-scale eddies (ICuzzi et al.u2008l) occur- 
ring near the dissipative scale of the turbulence. Larger m-sized 
particles pile up on large scales (i.e. larger than the gas scale 
height) in long-lived geost rophic pressure bumps surrounded by 
axisymmetr ic zonal flows (Dohansen et"aDj2009ab . In the model 
presented in lJohansen et aLl feOO?*) Thereafter referred to as J07], 
approximately meter- sized particles settle to form a thin mid- 
plane layer in balance between sedimentation and stirring by the 
gas which has de veloped turbulence throu gh the magnetorota- 
tional instability (iBalbus & Hawleyl Il99ll) . Particles then con- 
centrate in nearly axisymmetric gas high-pr essure regions which 
appear spontaneously in the turbulent flow (iFromang & Nelsonl 
120051: IJohansen et all 120061: iLvra et al.L |2008a|), reaching local 
column densities up to ten times the average. The passive con- 
centration is augmented as particles locally accelerate the gas 
towards the Keplerian speed, which leads to accumulation of 
particles drifting rapidly in from exterior orbits (a manifesta - 
tion of the streaming instability of lYoudin & Goodmani l2005l) . 
The gravitational attraction between the particles in the over- 
dense regions becomes high enough to initiate first a slow ra- 
dial contraction, and as the local mass density becomes com- 
parable to the Roche density, a full non-axisymmetric collapse 
to form gravitationally bound clumps with masses comparable 
to the 950-km-diameter dwarf planet Ceres (Mceres ~ 9.4 x 
10^^ kg). Such large pianetesimai birth sizes are in agreement 
with constraints from the current observed size distribution of 
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the asteroid belt (iMorbidelli et all l2Q09h and Neptune Trojans 
(ISheppard & TruiilloLl20T^ 



Some of the open questions related to this picture of plan- 
etesim al formation is to what degree the results of lJohansen et alJ 
(I2OO7I) are affected by the fact that self-gravity was turned on af- 
ter particles had concentrated in a pressure bumps and how the 
emergence and amplitude of pressure bumps are affected by nu- 
merical resolution. In this paper we present high-resolution and 
long-time-integration simulations of planetesimal formation in 
turbulence caused by the magnetorotational instability (MRI). 
We find that the large-scale geostrophic pressure bumps that are 
responsible for particle concentration are sustained when going 
from moderate (256^) to high (512^) resolution. Particle concen- 
tration in these pressure bumps is also relatively independent on 
resolution. We present a long-time-integration simulation per- 
formed at moderate resolution (256^) where particles and self- 
gravity are started at the same time, in contrast to earlier simu- 
lations where self-gravity was not turned on until a strong con- 
centration event occurred (J07). We also study the initial burst of 
planetesimal formation at 512^ resolution. We present evidence 
for collisions between gravitationally bound clumps, observed 
at both moderate and high resolution, and indications that the 
Initial Mass Function of gravitationally bound clumps involves 
masses ranging from a significant fraction of to several times 
the mass mass of Ceres. We point out that the physical nature 
of the collisions is unclear, since our numerical algorithm does 
not allow clumps to contract below the grid size. Gravitational 
scattering and binary formation are other possible outcomes of 
the close encounters, in case of resolved dynamics. Finding the 
Initial Mass Function of planetesimals forming from the gravita- 
tionally bound clumps will ultimately require an improved algo- 
rithm for the dynamics and interaction of bound clumps as well 
as the inclusion of particle shattering and coagulation during the 
gravitational contraction. 

The paper is organised as follows. In Sect. [2] we describe the 
dynamical equations for gas and particles. Sect. [3] contains de- 
scriptions of a number of improvements made to the Pencil Code 
in order to be able to perform particle-mesh simulations at up to 
at least 4096 cores. In Sect. |4] we explain the choice of simulation 
parameters. The evolution of gas turbulence and large-scale pres- 
sure bumps is analysed in Sect. [S] Particle concentration in sim- 
ulations with no self-gravity is described in Sect. [6l Simulations 
including particle self-gravity are presented in Sect. [7] (256^ res- 
olution) and Sect. [8] (5 12^ resolution). We summarise the paper 
and discuss the implications of our results in Sect. [9l 



2. Dynamical equations 

We perform simulations solving the standard shearing box 
MHD/drag force/self-gravity equations for gas defined on a fixed 
grid and solid particles evolved as numerical superparticles. We 
use the Pencil Code, a sixth order spatial and third order tempo- 
ral symmetric finite difference codeQ. 

We model the dynamics of a protoplanetary disc in the 
shearing box approximation. The coordinate frame rotates at 
the Keplerian frequency Q at an arbitrary distance ro from the 
central star. The axes are oriented such that the x points radi- 
ally away from the central gravity source, y points along the 
Keplerian flow, while z points vertically out of the plane. 



2.1. Gas velocity 

The equation of motion for the gas velocity u relative to the 
Keplerian flow is 

5m ^ „^ (0)du ^ ^ 1 ^ ^ ^ 

— -h (m • v)m -h w;^— = IQuyCx Uuxey -h zUAvex 

ot oy '2 

+ -Jx{B + Boz) - -VP - ^^(u - V) + /,(«) . (1) 

P p Tf 

The left hand side includes advection both by the velocity field 
u itself and by the linearised Keplerian flow uf^ = -(3/2)Qx. 
The first two terms on the right hand side represent the Coriolis 
force in the x- and ^-directions, modified in the j-component by 
the radial advection of the Keplerian flow, Uy = -u^ 
The third term mimics a global radial pressure gradient which 
reduces the orbital speed of the gas by the positive amount Av. 
The fourth and fifth terms in Eq. ([T]) are the Lorentz and pressure 
gradient forces. The current density is calculated from Ampere's 
law juqJ = V X B. The Lorentz force is modified to take into ac- 
count a mean vertical field component of strength ^o- The sixth 
term is a drag force term is described in Sect. 12 41 

The high-order numerical scheme of the Pencil Code has 
very little numerical diss i pation from time- stepping the advec- 
tion term (iBrandenburgi l2003l) . so we add explicit viscosity 
through the term fy(u) in Eq. ([T]). We use sixth order hypervis- 
cosity with a constant dynamic viscosity JI3 = vsp. 



/v 



Pg 



(2) 



This form of the viscosity conserves momentum. The opera- 
tor isJefinedasV^^^ -h d^/dy^ d^/dz^. It was shown 
by I Johansen et al.l (l2009ah that hyperviscosity simulations show 
zonal flows and pressure bumps very similar to simulations us- 
ing Navier-Stokes viscosity. 

2.2. Gas density 

The continuity equation for the gas density p is 



^ -h (m • V)p -h u"^^^ = -pV • u -h foip) . 



(0)^ 

dt ' ' ■ dy 



The diffusion term is defined as 
/d = Z)3VV, 



(3) 



(4) 



where is the hyperdiffusion coeflficient necessary to suppress 
Nyquist scale wiggles arising in regions where the spatial den- 
sity variation is high. We adopt an isothermal equation of state 
with pressure P = c^p and (constant) sound speed c^. 

2.3. Induction equation 

The induction equation for th e magnetic vector potential A is 
(see iBrandenburg et al.L 1 1 9951 for details) 



dA 

dt 



,/0) ' 

y dy 



+ wr— =ux(B-\- Boz) + -QAyX -h f^(A) . 



(5) 



See |http : //code . google . com/p/pencil-code7] 



The resistivity term is 

/, = ^3V^A, (6) 

where 773 is the hyperresistivity coeflScient. The magnetic field is 
calculated from B = V x A. 



Johansen et aL: High-resolution simulations of planetesimal formation 



3 



2.4. Particles and drag force scheme 



2.7. Collisions 



The dust component is treated as a number of individual super- 
particles, the position x and velocity v of each evolved according 
to 



djc 3 ^ 

-=v--Qxey., 

dv 1 



dt y X ^ X y V 



(7) 
(8) 



The particles feel no pressure or Lorentz forces, but are subjected 
to the gravitational potential of their combined mass. Particle 
collisions are taken into account as well (see Sect. l2.7l below). 

Two-way drag forces between gas defined on a fixed grid 
and Lagrangian particles are calculated through a particle-mesh 
method (see Youdin & Johansen, 2007, for details). First the gas 
velocity field is interpolated to the position of a particle, using 
second order spline interpolation. The drag force on the particle 
is then trivial to calculate. To ensure momentum conservation 
we then take the drag force and add it with the opposite sign 
among the 27 nearest grid points, using the Triangular Shaped 
Cloud scheme to ensure momentum conservation in the assign- 
ment (iHocknev & Eastwood Il98lh . 

2.5. Units 

All simulations are run with natural units, meaning that we set 
= Q = jdQ = Pq = I. Here po represents the mid-plane gas 
density, which in our unstratified simulations is the same as the 
mean density in the box. The time and velocity units are thus 
[t] = and [v] = c^. The derived unit of the length is the 
scale height H = c^jQ = [/]. The magnetic field unit is [B] = 

2.6. Self -gravity 

The gravitational attraction between the particles is calculated 
by first assigning the particle mass density pp on the grid, using 
the Triangular Shaped Cloud scheme described above. Then the 
gravitational potential at the grid points is found by inverting 
the Poisson equation 



V^0 = 47rGpp 



(9) 



using a Fast Fourier Transform (FFT) method (see online sup- 
plement of J07). Finally the self-potential of the particles is in- 
terpolated to the positions of the particles and the acceleration 
added to the particle equation of motion (Eq.[8]). We define the 
strength of self-gravity through the dimensionless parameter 



47tG 



(10) 



where po is the gas density in the mid-plane. This is practical 
since all simulations are run with ^2 = po = Cs = 1. Using 
i7g = ^/2npoH for the gas column density, we obtain a connec- 
tion between the dimensionless G and the relevant dimensional 
parameters of the box. 

We assume a standard ratio of particle to gas column densities of 
0.01. The self-gravity of the gas is ignored in both the Poisson 
equation and the gas equation of motion. 



Particle collisions become important inside dense particle 
clumps. In J07 the eff'ect of particle collisions was included in 
a rather crude way by reducing the relative rms speed of par- 
ticles inside a grid cell on th e collisiona l time- scale, to mimic 
coUisional cooling. Recently iRein et al" l (l20Toh found that the 
inclusion of particle collisions suppresses condensation of small 
scale clumps from the turbulent flow and favours the formation 
of larger structures. 

iRein et al.l (1201 Ol) claimed that the lack of collisions is 
an inherent flaw i n the superparticle approach. However, 
iLithwick & ChiangI (l2007b presented a scheme for the colli- 
sional evolution of particle rings whereby the collision between 
two particles occupying the same grid cell is determined by 
drawing a random number (to determine whether the two par- 
ticles are at the same vertical position). We have extended this 
algorithm to model collisions between superparticles based on 
a Monte Carlo algorithm. We obtain the correct collision fre- 
quency by letting nearby superparticle pairs collide on the av- 
erage once per collisional time-scale of the swarms of physical 
particles represented by each superparticle. 

We have implemented the superparticle collision algorithm 
in the Pencil Code and will present detailed numerical tests that 
show its validity in a paper in preparation (Johansen, Youdin, 
& Lithwick, in preparation). The algorithm gives each particle a 
chance to interact with all other particles in the same grid cell. 
The characteristic time-scale r^^^^^ = l/(njCrijSvij) for a represen- 
tative particle in the swarm of superparticle / to collide with any 
particle from the swarm represented by superparticle j is calcu- 
lated by considering the number density nj represented by each 
superparticle, the collisional cross section ct/j of two swarm par- 
ticles, and the relative speed Svij of the two superparticles. For 
each possible collision a random number P is chosen. If P is 
smaller than St/Zcoiu where St is the time-step set by magne- 
tohydrodynamics, then the two particles collide. The collision 
outcome is determined as if the two superparticles were actual 
particles with radii large enough to touch each other. By solving 
for momentum conservation and energy conservation, with the 
possibility for inelastic collisions to dissipate kinetic energy to 
heat and deformation, the two colliding particles acquire their 
new velocity vectors instantaneously. 

All simulations i nclude c ollisions with a coefficient of resti- 
tution of 6 = 0.3 (e.g lBlum"& Muench, 1992), meaning that each 
collision leads to the dissipation of approximately 90% of the 
relative kinetic energy to deformation and heating of the collid- 
ing boulders. We include particle collisions here to obtain a more 
complete physical modelling. Detailed tests and analysis of the 
effect of particle collisions on clumping and gravitational col- 
lapse will appear in a future publication (Johansen, Youdin, & 
Lithwick, in preparation). 



3. Computing resources and code optimisation 

For this project we were kindly granted access to 4096 cores on 
the "Jugene" Blue Gene/P system at the Jiilich Supercomputing 
Centre (JSC) for a total of five months. Each core at the 
BlueGene/P has a clock speed of around 800 MHz. The use of 
the Pencil Code on several thousand cores required both trivial 
and more fundamental changes to the code. We describe these 
technical improvements in detail in this appendix. 

In the following nxgrid, nygrid, nzgrid refer to the full 
grid dimension of the problem. We denote the processor num- 
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512^ gas + 64x10^ particles 
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Fig. 1. Scaling test for particle-mesh problem with 512^ grid 
cells and 64 x 10^ particles. The particles are distributed evenly 
over the grid, so that each core has the same number of particles. 
The inverse code speed is normalised by the number of time- 
steps and by either the total number of grid points and particles 
(top panel) or by the number of grid points and particles per core 
(bottom panel). 
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Fig. 2. Code speed as a function of simulation time (top panel) 
and maximum particle number on any core (bottom panel) for 
256^ resolution on 2048 cores. Standard domain decomposition 
(SDD) quickly becomes unbalanced with particles and achieves 
only the speed of the particle-laden mid-plane cores. With the 
Particle Block Domain Decomposition (PBDD) scheme the 
speed stays close to its optimal value, and the particle number 
per core (bottom panel) does not rise much beyond 10^. 



ber by ncpus and the directional processor numbers as nprocx, 
nprocy, nprocz. 

3.1. Changes made to the Pencil Code 
3.1.1. Memory optimisation 

We had to remove several uses of global arrays, i.e. 2-D or 
3-D arrays of linear size equal to the full grid. This mostly 
affected certain special initial conditions and boundary condi- 
tions. An additional problem was the use of an array of size 
(ncpus , ncpus) in the particle communication. This array was 
replaced by a 1-D array with no problems. 

The runtime calculation of 2-D averages (e.g. gas and par- 
ticle column density) was done in such a way that the whole 
(nxgrid,nygrid) array was collected at the root processor in 
an array of size (nxgrid , nygrid , ncpus) , before appending a 
chunk of size (nxgrid, nygrid) to an output file. The storage 
array, used for programming convenience in collecting the 2-D 
average from all the relevant cores, became excessively large at 
high resolution and processor numbers, and we abandoned the 
previous method in favour of saving chunks of the column den- 
sity array into separate files, each maintained by the root proces- 



sors in the z-direction. A similar method was implemented for 
}^-averages and x- averages. 

The above-mentioned global arrays had been used in the 
code for programming convenience and did not imply excessive 
memory usage at moderate resolution and processor numbers. 
Purging those arrays in favour of loops or smaller 1-D arrays 
was relatively straight-forward. 

3.1.2. Particle migration 

At the end of a sub-time-step each processor checks if any par- 
ticles have left its spatial domain. Information about the number 
of migrating particles, and the processors that they will migrate 
into, is collected at each processor. The Pencil Code would then 
let all processors exchange migrating particles with all other pro- 
cessors. In practice particles would of course only migrate to 
neighbouring processors. However, at processor numbers of 512 
or higher, the communication load associated with each proces- 
sor telling all other processors how many particles it wanted to 
send (mostly zero) was so high that it dominated over both the 
MHD and the particle evolution equations. Since particles in 
practice only migrate to the neighbouring processors any way. 
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we solved this problem by letting the processors only communi- 
cate the number of migrating particles to the immediate neigh- 
bours. Shear-periodic boundary conditions require a (simple) al- 
gorithm to determine the three neighbouring processors over the 
shearing boundary in the beginning of each sub-time- step. 

3.2. Timings 



With the changes described in Sect. 13.1.11 and Sect. 13.1.21 the 
Pencil Code can be run with gas and particles efficiently at sev- 
eral thousand processors, provided that the particles are well- 
mixed with the gas. 

In Fig. [T] we show timings for a test problem with 512^ grid 
cells and 64 x 10^ particles. The particles are distributed evenly 
over the grid, avoiding load balancing challenges described be- 
low. We evolve gas and particles for 1000 time-steps, the gas 
and particles subject to standard shearing box hydrodynamics 
and two-way drag forces. The lines show various drag force 
schemes - NGP correspondi ng to Nearest Grid Point and TSC 
to Triangular Shaped Cloud (iHockney & Eastwood 1 198 ll) . We 
achieve near perfect scaling up to 4096 cores. Including self- 
gravity by a Fast Fourier Transform (FFT) method the code 
slows down by approximately a factor of two, but the scaling 
is only marginally worse than optimal, with less than 50% slow- 
down at 4096 cores. This must be seen in the light of the fact that 
3-D FFTs involve several transpositions of the global density ar- 
ray, each transposition requiring the majority of grid points to be 
communicated to another processor (see online supplement of 
J07). 

3.3. Particle parallelisation 

At high resolution it becomes increasingly more important to 
parallelise efficiently along at least two directions. In earlier pub- 
lications we had run 256^ simulations at 64 processors, with the 
domain decomposition nprocy=32 and nprocz=2 (J07). Using 
two processors along the z-direction exploits the intrinsic mid- 
plane symmetry of the particle distribution, while the Keplerian 
shear suppresses most particle density variation in the azimuthal 
direction, so that any processor has approximately the same 
number of particles. 

However, at higher resolution we need to either have 
nprocz>2 or nprocx>l, both of which is subject to particle 
clumping (from either sedimentation or from radial concentra- 
tions). This would in some cases slow down the code by a factor 
of 8-10. We therefore developed an improved particle paralleli- 
sation, which we denote Particle Block Domain Decomposition 
(PBDD). This new algorithm is described in detail in the follow- 
ing subsection. 

3.3.1. Particle Block Domain Decomposition 

The steps in Particle Block Domain Decomposition scheme are 
as follows: 

1 . The fixed mesh points are domain-decomposed in the usual 
way (with ncpus=nprocxxnprocyxnprocz). 

2. Particles on each processor are counted in bricks of size 
nbxxnbyxnbz (typically nbx=nby=nbz=4). 

3. Bricks are distributed among the processors so that each pro- 
cessor has approximately the same number of particles 

4. Adopted bricks are referred to as blocks. 

5. The Pencil Code uses a third order Runge-Kutta time- 
stepping scheme. In the beginning of each sub-time- step par- 



ticles are counted in blocks and the block counts commu- 
nicated to the bricks on the parent processors. The particle 
density assigned to ghost cells is folded across the grid, and 
the final particle density (defined on the bricks) is commu- 
nicated back to the adopted blocks. This step is necessary 
because the drag force time- step depends on the particle den- 
sity, and each particle assigns density not just to the nearest 
grid point, but also to the neighbouring grid points. 

6. In the beginning of each sub-time- step the gas density and 
gas velocity field is communicated from the main grid to the 
adopted particle blocks. 

7. Drag forces are added to particles and back to the gas grid 
points in the adopted blocks. This partition aims at load bal- 
ancing the calculation of drag forces. 

8. At the end of each sub-time-step the drag force contribution 
to the gas velocity field is communicated from the adopted 
blocks back to the main grid. 

To illustrate the advantage of this scheme we show in Fig. [2] 
the instantaneous code speed for a problem where the particles 
have sedimented to the mid-plane of the disc. The grid resolution 
is 256^ and we run on 2048 cores, with 64 cores in the j-direction 
32 cores in the z-direction. The blue (black) line shows the re- 
sults of running with standard domain decomposition, while the 
orange (grey) line shows the speed with the improved Particle 
Block Domain Decomposition scheme. Due to the concentration 
of particles in the mid-plane the standard domain decomposition 
leaves many cores with few or no particles, giving poor load 
balancing. This problem is alleviated by the use of the PBDD 
scheme (orange/grey line). 

PBDD works well as long as the single blocks do not achieve 
higher particle density than the optimally distributed particle 
number npar/ncpus. In the case of strong clumping, e.g. due 
to self-gravity, the scheme is no longer as efficient. In such ex- 
treme cases one should ideally limit the local particle number in 
clumps by using sink particles. 



4. Simulation parameters 

The main simulations of the paper focus on the dynamics and 
self-gravity of solid particles moving in a gas flow which has 
developed turbulence thro ugh the magnetorotational instability 
(iBalbus & Hawleyl Il99ll) . We have performed two moderate- 
resolution simulations (with 256^ grid points and 8 x 10^ par- 
ticles) and one high-resolution simulation (512^ grid points and 
64 X 10^ particles). Simulation parameters are given in Table [T] 
We use a cubic box of dimensions (1.32//)^. 

Note that we use a sub-Keplerian speed difl'erence of Av = 
0.05 which is higher than Av = 0.02 presented in the main paper 
of J07. The ability of pressure bumps to trap particles is gener- 
ally reduced with increasing Av (see online supplementary in- 
formation for J07). Particle clumping by streamin g instabilities 
also becomes less efficient as Av is increased (JOT. lBai & Stonel 
2010c). Estimates of the radial pressure support in discs can 
be extracted from models of column density and temperature 
structure. A gas parcel orbiting at a radial distance r from the 
star, where the disc aspect ratio is H/r and the mid-plane radial 
pressure gradient is d In P/d In r, orbits at a sub-Keplerian speed 
V = vk - Av. The speed reduction Av is given by 



Av 



IHdlnP 
2 r dlnr 



(12) 



In the Minimum Mass Solar Nebula of 'Havashi|__( 19811) 



din P/d In r = -3.25 in the mid-plane (e.g. lYoudia,2008h . The 
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Table 1. Simulation parameters in natural units. 



Run 


L^xLyX 


N^xNyX Nz 


V3 = Ds 


= V3 


Bo 




Av 


G 


A^ 


^par 


^grav,l 


^grav,2 


Ml 


(1.32)3 


2563 


2x 10- 


-14 


0.0015 


9x 10^ 


0.05 


N/A 


40.0 


N/A 


N/A 


N/A 


M2 


(1.32)3 


2563 


2x 10- 


-14 


0.003 


2x 10^ 


0.05 


0.5 


40.0 


20.0 


20.0 


33.0 


H 


(1.32)3 


5123 


7x 10- 


-16 


0.0015 


9x 10^ 


0.05 


0.5 


40.0 


32.0 


35.0 


37.0 



The first column gives the name of the simulation, while the box size and the grid resolution are given in the following two columns. The next 
column gives the hyperdiff'usivity coefficients. The next two columns give the mean vertical magnetic field and the corresponding plasma-/?. The 
next column gives the sub-Keplerian velocity difference. The following column shows the particle self-gravity parameter G. The last four columns 
give the total simulation time, the time when particles were released, and the times at which self-gravity was started and self-gravity simulations 
were stopped. 



resulting scaling of the sub-Keplerian speed with the orbital dis- 
tance is 



Av 



/ r \ 



1/4 



(13) 



The slightly colder disc model used by iBrauer et"ar ] (I2008ah 
yields instead 



— ^0.04 — 

Cc \AU/ 



(14) 



In more complex disc models the pressure profile is changed 
e.g. in in terfaces between regions of weak and strong turbulence 
(iLyra et a L, 200 8^. We use throughout this paper a fixed value 
of Av/cs = 0.05. 

Another diflTerence from the simulations of JOT is that the 
turbulent viscosity of the gas flow is around 2-3 times higher, 
because of the increased turbulent viscosity when going from 
256^ to 512^ (see Sect. O. Therefore we had to use a stronger 
gravity in this paper, G = 0.5 compared to G = 0.2 in JOT, to see 
bound particle clumps (planetesimals) forming. We discuss the 
implications of using a higher disc mass further in our conclu- 
sions. 

In all simulations we divide the particle component into four 
size bins, with friction time Qzf = 0.25, 0.50, 0.T5, 1.00, respec- 
tively. The particles drift radially because of the headwind from 
the g as orbiting at velocity Uy = -Av relative to the Keplerian 
flow (IWeidenschillingl [l9TTah . It is important to consider a dis- 
tribution of particle sizes, since the dependence of the radial drift 
on the particle sizes can have a negative impact on the ability 
of the particle mid-plane layer to undergo gravitational collapse 
dWeidenschilling, 1995). 

The Minimum Mass Solar Nebula has column den sity i7g = 
lT00(r/AU)-i-5gcm-2 ^ 150 gcm'^ at r = 5AU (iHavashiL 
Il98lh . and thus G = 0.042 according to Eq. O- Since we 
use ~10 times higher value for G, the mean free path of gas 
molecules is only /l~10cm. Therefore our choice of dimen- 
sionless friction times Qzf = 0.25, 0.50, 0.T5, 1.00 puts parti - 
cles in the Stokes drag force regime (IWeidenschillingl ll9TTal) . 
Here the friction time is independent of gas density, and the 
Stokes number Qzf is proportional to particle radius squared, 
so Qzf = 0.25, 0.50, 0.T5, 1.00 correspond to physical particle 
sizes ranging from 40 cm to 80 cm (see online supplement of 
JOT). Scaling Eq. (fTTI) to more distant orbital locations gives 
smaller physical particles and a gas column density closer to the 
Minimum Mass Solar Nebula value, since self-gravity is more 
efiftcient in regions where the rotational support is lower. 

There are several points to be raised about our protoplanetary 
disc model. The high self-gravity parameter that we use implies 
not only a very high column density, but also that the gas compo- 
nent is close to gravitational instability. The self-gravity parame- 
ter G is connected to the more commonly used Q (where 2 > 1 is 



the axisyr ametric stability criterion for a fla t disc in Keplerian ro- 
tation, see ISafronovl 1 1 960l: iToomrel 1 1 964l) through G ~ 1.62"^. 
Thus we have 2 3.2, which means that gas self-gravity should 
start to aflTect the dynamics (the disc is not formally gravitation- 
ally unstable, but the disc is massive enough to slow down the 
propagation of sound waves). Another issue with such a mas- 
sive disc is our assumption of ideal MHD. The high gas col- 
umn density decreases the ionisation by cosmic rays and X-rays 
and increases the recombinati on rate on dust grains (ISano et al.L 
120001: iFromang et aP. l2002h . iLesur & Longarettil (l200Tl l201ol) 
have furthermore shown that the ratio of viscosity to resistiv- 
ity, the magnetic Prandtl number, afiTects both small-scale and 
large-scale dynamics of saturated magnetorotational turbulence. 
Ideally all these eflTects should be taken into account. However, 
in this paper we choose to focus on the dynamics of solid par- 
ticles in gas turbulence. Thus we include many physical eflTects 
that are important for particles (drag forces, self-gravity, colli- 
sions), while we ignore many other eflTects that would determine 
the occurrence and strength of gas turbulence. This approach al- 
lows us to perform numerical experiments which yield insight 
into planetesimal formation with relatively few free parameters. 

4.1. Initial conditions 



The gas is initialised to have unit density everywhere in the box. 
The magnetic field is constant B = Boe^. The gas velocity field 
is set to be sub-Keplerian with Uy = -Av, and we furthermore 
perturb all components of the velocity field by small random 
fluctuations with amplitude Sv = 0.001, to seed modes that are 
unstable to the magnetorotational instability. In simulations with 
particles we give particles random initial positions and zero ve- 
locity. 



5. Gas evolution 

We start by describing the evolution of gas without particles, 
since the large-scale geostrophic pressure bumps appearing in 
the gas controls particle concentration and thus the overall abil- 
ity for planetesimals to form by self-gravity. The most important 
agent for dri ving gas dynamics is the magnetorotational insta- 
bility (MRI, iBalbus & HawlevL Il99lh which exhibits dynami- 
cal growth when the vertical component of the magnetic field 
is not too weak or too strong. The non-stratified MRI saturates 
to a s tate of non-linear subsonic fluctuations (e.g.'Ha wley et all 
fT995h . In this state there is an outward angular momentum flux 
through hydrodynamical Reynolds stresses (pUxUy) and magne- 
tohydrodynamical Maxwell stresses {-[i^^B^By). 

In Fig. [3] we show the Maxwell and Reynolds stresses as a 
function of time. Using a mean vertical field of = 0.0015 
(corresponding to plasma-beta oi p ^ 9 x 10^) the turbu- 
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-0.6-0.4-0.2 0.0 0.2 0.4 0.6 -0.6-0.4-0.2 0.0 0.2 0.4 0.6 -0.6-0.4-0.2 0.0 0.2 0.4 0.6 
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Fig. 4. The gas column density averaged over the azimuthal direction, as a function of radial coordinate x and time t in orbits. 
Large-scale pressure bumps appear with approximately 1% amplitude at both 256^ and 512^ resolution. 



lent viscosity almost triples when going from 256^ to 512^ 
grid points. This is in stark contrast with zero net flux simula- 
tions that show decreasing turb ulence with increasing resolution 
(iFromang & Papaloizoul 120071) . We interpret the behaviour of 
our simulations as an efl'ect of underresolving the most unstable 
wavelength of the magnetorotational instability. Considering a 
vertical magnetic field of const ant strength ^p, the most unsta- 
ble wave number of the MRI is (lBalbus&HawlevLll99lh 



(15) 




where va = BqI ^JiiqPq is the Alfven speed. The most un- 
stable wavelength is = Injkz. For = 0.0015 we get 
Az ~ 0.01//. The resolution elements are 6x ^ 0.005// at 256^ 
and Sx ^ 0.0026// at 512^. Thus we get a significant improve- 
ment in the resolution of the most unstable wavelen gth when go- 
ing fr om 256^ to 512^ g rid points. Other authors (ISimon et al.L 
[20091: F^i^ et al.L [2Q09f) have reported a similar increase in tur- 
bulent activity of net-flux simulations with increasing resolution. 
Our simulations show that this increase persists up to at least 
yS^9xlO^ 

The original choice of = 0.0015 was made in J07 in 
order to prevent the turbulent viscosity from dropping below 
a = 0.001. However, we can not obtain the same turbulent vis- 
cosity (i.e. a ~ 0.001) at higher resolution, given the same ^o- 
For this reason we did all 256^ experiments on particle dynamics 
and self-gravity with = 0.003 (J3 ^ 2x10^), yielding approx- 



imately the same turbulent viscosity as in the high-resolution 
simulation. 

The Reynolds and Maxwell stresses can be translated 
into an average turbulent viscosity (following the notation of 
Brandenburg et al., 1995), 



2 



1 



(16) 
(17) 



Here v, 



-^^"^ and 



,(mag) 



are the turbulent viscosities due to the ve- 



locity field and the magnetic field, respectively. We can further 
normalise the turbulent viscosities by the so und speed and gas 
scale height H (Shakur a & Sunyaev.i 19731) . 



1 r (kin) , (mag)1 



(18) 



0.0022, and 



We thus find a turbulent viscosity ofa ^ 0.001, a 
a ^ 0.003 for runs Ml, M2, and H, respectively. 

The combination of radial pressure support and two- 
way drag forces allows systematic relative motion be- 
tween gas an d particles, which is unstable to the stream- 



ing instabilitv (lYoudin & Goodniani 12005 : lYoudir 
'2007; "Jo hansen & YoudinL l2007l: iMiniatl l2oTQr 



& Johansen. 



Bai & StoneL 



JOlOa^b). Streaming instabilities and ma gnetorotationa l insta- 
bilities can opera te in concurrence (J07 iBalsara et al.L l2009l: 
iTillev etal.L[20TQh . However, we find that particles concentrate 
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Fig. 3. Maxwell and Reynolds stresses as a function of time. 
The Reynolds stress is approximately five times lower than the 
Maxwell stress. There is a marked increase in the turbulent 
stresses when increasing the resolution from 256^ to 512^ at a 
fixed mean vertical field = 0.0015, likely due to better reso- 
lution of the most unstable MRI wavelength at higher resolution. 
Using ^0 = 0.003 at 256^ gives turbulence properties more sim- 
ilar to 5 12l 



in high-pressure bumps forming due to the MRI, so that stream- 
ing instabilities are a secondary eff'ect in the simulations. A ne- 
cessity for the streaming instability to operate is a solids-to-gas 
ratio that is locally at least unity. The particle density in the mid- 
plane layer is reduced by turbulent diff'usion (which is mostly 
caused by the MRI), so in this way an increase in the strength 
of MRI turbulence can reduce the importance of the SI. Even 
though streaming instabilities do not appear to be the main driver 
of particle concentration in our simulations, the back-reaction 
drag force of the particles on the gas can potentially play a role 
during the gravitational contraction phase where local particle 
column densities get very high. The high gas column density 
needed for gravitational collapse in the current paper may also 
in reality preclude activity by the magnetorotational instability, 
given the low ionisation level in the mid-plane, which would 
make the streaming instability the more likely path to clumping 
and planetesimal formation. 

5.1. Pressure bumps 

An important feature of magnetorotational turbulence is the 
emergence of large-scale slowly overturning p ressure bumps 
(iFromang & NelsonLl2005l:ljohansen et aLLl2006l) . Such pressure 
bumps form with a zonal flow envelope due to random excitation 
of the zonal flow by lar ge-scale variations in the Maxwell stress 
(I Johansen et aLL[2009ah . Variations in the mean field magnitude 
and direction has also been shown to lead to the formation of 
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Fig. 5. The particle column density averaged over the azimuthal 
direction, as a function of radial coordinate x and time t in orbits. 
The starting time was chosen to be slightly prior to the emer- 
gence of a pressure bump (compare with left-most and right- 
most plots of Fig. [4]). The particles concentrate slightly down- 
stream of the gas pressure bump, with a maximum column den- 
sity between three and four times the mean particle column den- 
sity. The particles are between 40 and 80 cm in radius (i.e. boul- 
ders) for our adopted disc model. 



pressure bu mps in the interface region between weak and strong 
turbulence (iKato et al.|j 2009. 2010>) . Pressure bumps can also be 
launched b y a radial variation in resistivity, e.g. at the e dges of 
dead zones (iLvra et al.L l2008bl: iDzvurkevich et alll201Qh . 

Large particles - pebbles, rocks, and boulders - are attracted 
to the center of pressure bumps, because of the drag force asso- 
ciated with the sub-Keplerian/super-Keplerian zonal flow enve- 
lope. In presence of a mean radial pressure gradient the trapping 
zone is slightly downstream of the pressure bump, where there 
is a local maximum in the combined pressure. 

An efficient way to detect pressure bumps is to average the 
gas density field over the azimuthal and vertical directions. In 
Fig. |4] we show the gas column density in the 256^ and 512^ 
simulations averaged over the j^-direction, as a function of time. 
Large-scale pressure bumps are clearly visible, with spatial cor- 
relation times of approximately 10-20 orbits. The pressure bump 
amplitude is around 1%, independent of both resolution and 
strength of the external field. Larger boxes have been shown 
to result in higher- amp litude and longer-lived pressure bumps 
(I Johansen et al.Ll2009ah . We limit ourselves in this paper to a rel- 
atively small box, where we can achieve high resolution of the 
gravitational collapse, but plan to model planetesimal formation 
in larger boxes in the future. 
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Fig. 6. The particle column density as a function of time after self-gravity is turned on at ^ = 20.0rorb, for run M2 (256^ grid cells 
with 8x10^ particles). Each gravitationally bound clump is labelled by its Hill mass in units of Ceres masses. The insert shows an 
enlargement of the region around the most massive bound clump. The most massive clump at the end of the simulation contains a 
total particle mass of 34.9 Ceres masses, partially as the result of a collision between a 13.0 and a 14.6 Ceres mass clump occurring 
at a time around t = 31.6rorb- 
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Fig. 7. Temporal zoom in on the collision between three clumps (planetesimals) in the moderate resolution run M2. Two clumps 
with a radial separation approximately 0.05H shear past each other, bringing their Hill spheres in contact (first two panels). The 
clumps first pass through each other (panels three and four), but eventually merge (fifth panel). Finally a much lighter clump collides 
with the newly formed merger product (sixth panel). 
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6. Particle evolution 

We release the particles at a time when the turbulence has satu- 
rated, but choose a time when there is no significant large-scale 
pressure bump present. Thus we choose t = 20rorb for the 256^ 
simulation and t = 32rorb for the 512^ simulation (see left-most 
and right-most plot of Fig. [4]). In the particle simulations we al- 
ways use a mean vertical field Bq = 0.003 at 256^ to get a turbu- 
lent viscosity more similar to 512^. The four friction time bins 
(QTf = 0.25,0.50,0.75,1.00)correspondto particle sizes between 
40 and 80 cm. 

The particles immediately fall towards the mid-plane of the 
disc, before finding a balance between sedimentation and tur- 
bulent stirring. Fig. [5] shows how the presence of gas pressure 
bumps has a dramatic influence on particle dynamics. The par- 
ticles display column density concentrations of up to 4 times 
the average density just downstream of the pressure bumps. At 
this point the gas moves close to Keplerian, because the (posi- 
tive) pressure gradient of the bump balances the (negative) ra- 
dial pressure gradient there. The column density concentration 
is relatively independent of the resolution, as expected since the 
pressure bump amplitude is almost the same. 



7, Self-gravity - moderate resolution 

In the moderate-resolution simulation (256^) we release parti- 
cles and start self-gravity simultaneously at ^ = 20rorb. This is 
difl'erent from the approach taken in J07 where self-gravity was 
turned on after the particles had concentrated in a pressure bump. 
Thus we address concerns that the continuous self-gravity inter- 
action of the particles would stir up the particle component and 
prevent the gravitational collapse. After releasing particles we 
continue the simulation for another thirteen orbits. Some repre- 
sentative particle column density snapshots are shown in Fig.[6j 
As time progresses the particle column density increases in high- 
pressure structures with power on length scales ranging from a 
few percent of a scale height to the full radial domain size. Self- 
gravity becomes important in these overdense regions, so some 
local regions begin to contract radially under their own weight, 
eventually reaching the Roche density and commencing a fully 
2-D collapse into discrete clumps. 

The Hill sphere of each bound clump is indicated in Fig. [6l 
together with the mass of particles encompassed inside the Hill 
radius (in units of the mass of the dwarf planet Ceres). We cal- 
culate the Hill radius of clumps at a given time by searching for 
the point of highest particle column density in the x-y plane. We 
first consider a "circle" of radius one grid point and calculate the 
two terms relevant for determination of the Hill radius - the tidal 
term 3Q^R and the gravitational acceleration GM^/R^ of a test 
particle at the edge of the Hill sphere due to the combined grav- 
ity of particles inside the Hill sphere. The mass Mp contained in 
a cylinder of radius R must fulfil 



3Q^R^ 



P G 



(19) 



The natural constant G is set by the non-dimensional form of the 
Poisson equation. 



j2 ^ _ 47rG Pp 



HereV^ = /d(x/H)^+d^ /d(y/H)^+d^ /d(z/H)^. Using nsitmsil 
units for the simulation, with = Q = H = po = I, together 
with our choice of 



47rG 



(21) 



we obtain an expression for the gravitational constant G. We then 
check the validity of the expression 



M^ = Y,n,PtjSV> 



UnpoR^ 



(22) 



where p^y = YukPijk is the vertically averaged mass density at 
grid point (/, j) and dV is the volume of a grid cell. It is conve- 
nient to work with p^^ since this vertical average has been output 
regularly by the code during the simulation. The sum in Eq. (l22l) 
is taken over all grid points lying within the circle of radius R 
centred on the densest point. We continue to increase R in units 
of dx until the boundness criterion is no longer fulfilled. This de- 
fines the Hill radius R. Strictly speaking our method for finding 
the Hill radius is only valid if the particles were distributed in a 
spherically symmetric way. In reality particles are spread across 
the mid-plane with a scale height of approximately 0.04//. We 
nevertheless find by inspection that the calculated Hill radii cor- 
respond well to the regions where the particle flow appears dom- 
inated by the self-gravity rather than the Keplerian shear of the 
main flow and that the mass within the Hill sphere does not fluc- 
tuate because of the inclusion of non-bound particles. 

The particle-mesh Poisson solver based on Fast Fourier 
Transforms can not consider the gravitational potential due to 
structure within a grid cell. From the perspective of the gravity 
solver the smallest radius of a bound object is thus the grid cell 
length 6x. This puts a lower limit to the mass of bound structures, 
since the Hill radius can not be smaller than a grid cell. 



^ /GM^y/3 
This gives a minimum mass of 



(23) 



(24) 



(20) 



Using = 2.0 x 10^^ g, Hjr = 0.05 and 6x = 0.0052/^ 
(256^) or dx = 0.0026 (512^), we get the minimum mass of 
Minin ~ O.llMceres at 256^ and M^,^ ^ 0.014Mceres at 5121 
Less massive objects will inevitably be sheared out due to the 
gravity softening. 

Fig. [6] shows that a number of discrete particle clumps con- 
dense out of the turbulent flow during the 13 orbits that are run 
with self-gravity. The most massive clump has the mass of ap- 
proximately 35 Ceres masses at the end of the simulation, while 
four smaller clumps have masses between 2.4 and 4.6 Ceres 
masses. The smallest clumps are more than ten times more mas- 
sive than the minimum resolvable mass. 



7.1. Planetesimal collision 

The 35 Ceres masses particle clump visible in the last panel 
of Fig. [6] is partially the result of a collision between a 13.0 
and a 14.6 Ceres mass clump at a time around t - 31.6rorb- 
The collision is shown in greater detail in Fig. [71 The merg- 
ing starts when two clumps with radial separation of approxi- 
mately 0.05// shear past each other, bringing their Hill spheres 




Fig. 8. The particle column density as a function of time after self-gravity is turned on after t = 35.0rorb in the high-resolution 
simulation (run H with 512^ grid cells and 64 x 10^ particles). Two clumps condense out within each other's Hill spheres and 
quickly merge. At the end of the simulation bound clumps have masses between 0.5 and 7.5 Mceres- 
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Fig. 9. The total bound mass and the mass in the most massive gravitationally bound clump as a function of time. The left panel 
shows the result of the moderate-resolution simulation (M2). Around a time of SOJorb there is a condensation event that transfers 
many particles to bound clumps. Two orbits later, at 32rorb, the two most massive clumps collide and merge. The right panel shows 
the high-resolution simulation (H). The total amount of condensed material is comparable to M2, but the mass of the most massive 
clump is smaller. This may be a result either of increased resolution in the self-gravity solver or of the limited time span of the high- 
resolution simulation. The total particle mass for both resolutions is Mtot ~ 460 Mceres- Only around 10% of the mass is converted 
into planetesimals during the time- span of the simulations. 



256' 



512' 




Fig. 10. Histogram of clump masses after first production of bound clumps and at the end of the simulation. At moderate resolution 
(left panel) only a single clump condenses out initially, but seven orbits later there are five clumps, including the 30-h Mceres object 
formed by merging. At high resolution (right panel) the initial planetesimal burst leads to the formation of many sub-Ceres-mass 
clumps. The most massive clump is similar to what forms initially in the moderate-resolution run. 



in contact. The less massive clumps passes first directly through 
the massive clump, appearing distorted on the other side, before 
merging completely. A third clump collides with the collision 
product shortly afterwards, adding another 3.5 Ceres masses to 
the clump. 

The particle-mesh self-gravity solver does not resolve par- 
ticle self-gravity on scales smaller than a grid cell. The bound 
particle clumps therefore stop their contraction when the size 
of a grid cell is reached. This exaggerated size increases the 
collisional cross section of planetesimal artificially. The clumps 
behave aerodynamically like a collection of dm- sized particles, 
while a single dwarf planet sized would have a much longer fric- 



tion time. Therefore the planetesimal collisions that we observe 
are not conclusive evidence of a collisionally evolving size dis- 
tribution. Future convergence tests at extremely high resolution 
(1024^ or higher), or mesh refinement around the clumps, will 
be needed to test the validity of the planetesimal mergers. 

The system is however not completely dominated by the 
discrete gravity sources. A significant population of free parti- 
cles are present even after several bound clumps have formed. 
Those free particles can act like dynamical friction and thus 
allow close encounters to lead to merging or binary formation 
(iGoldreich et al.L l2002h . In the high-resolution simulation pre- 
sented in the next section we find clear evidence of a trailing 
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spiral density structure that is involved in the collision between 
two planetesimals. 

8. Self-gravity - high resolution 

In Section[6]we showed that particle concentration is maintained 
when going from 256^ to 512-^ grid cells. In this section we show 
that the inclusion of self-gravity at high resolution leads to rapid 
formation of bound clumps similar to what we observe at mod- 
erate resolution. Given the relatively high computational cost of 
self-gravity simulations we start the self-gravity at ^ = SSTorb 
in the 512^ simulation, three orbits after inserting the particles. 
The evolution of particle column density is shown in Fig.[8l Due 
to the smaller grid size bound particle clumps appear visually 
smaller than in the 256^ simulation. The increased compactness 
of the planetesimals can potentially decrease the probability for 
planetesimal collisions (Sect. 17.1b , which makes it important to 
do convergence tests. 

The high-resolution simulation proceeds much as the 
moderate-resolution simulation. Panels 1 and 2 of Fig. [8] show 
how overdense bands of particles contract radially under their 
own gravity. The increased resolution of the self-gravity solver 
allows for a number of smaller planetesimals to condense out as 
the bands reach the local Roche density at smaller radial scales 
(panel 3). Two of the clumps are born within each other's Hill 
spheres. They merge shortly after into a single clump (panel 5). 
This clump has grown to 7.5 Mceres at the end of the simulation, 
which is the most massive clump in a population of clumps with 
masses between 0.5 and 7.5 Mceres- 

Although we do not reach the same time span as in the low 
resolution simulation, we do observe two bound clumps col- 
liding. However, the situation is different, since the colliding 
clumps form very close to each other and merge almost imme- 
diately. An interesting aspect is the presence of a particle den- 
sity structure trailing the less massive clump. The gravitational 
torque from this structure can play an important role in the colli- 
sion between the two clumps, since the clumps initially appear to 
orbit each other progradely . This confirms the trend observed in 
I Johansen & Lacerd"al(l201Ql) for particles to be accreted with pro- 
grade angular momentum in the presence of drag forces, which 
can explain why the largest asteroids tend to spin in the pro- 
grade directioifl. The gravitational torque from the trailing den- 
sity structure would be negative in that case and cause the rela- 
tive planetesimal orbit to shrink. 

8.1. Clump masses 

In Fig. [9] we show the total mass in bound clumps as a function of 
time. Finding the physical mass of a clump requires knowledge 
of the scale height H of the gas, as that sets the physical unit 
of length. The self-gravity solver in itself only needs to know 
G = 47iG/(Q^/po), which is eff'ectively a combination of density 
and length scale. When quoting the physical mass we assume 
orbital distance r = 5 AU and aspect ratio H/r = 0.05. The total 
mass in particles in the box is Mtot = O.OH^g^sLxLy ~ 460 Mceres, 
with G = 0.5 and iJgas = 1800 g cm'^ from Eq. 0- 

In both simulations approximately 50 Mceres of particles are 
present in bound clumps at the end of the simulation. However, 

^ Prograde rotation is not readily acquired in standard numerical 
models where planetesimals accumulate in a gas-free environment, al- 
though planetary birth in rotating self-gravitating gas blobs was recently 
been put forward to explain the prograde rotation of planets dNavakshinl . 



self-gravity was turned on and sustained for very diff'erent times 
in the two diff'erent simulations. In the moderate-resolution sim- 
ulation most mass is bound in a single clump at the end of the 
simulation. The merger event discussed in Sect. 17.11 is clearly 
visible around t = 32rorb. 

Fig. [TOl shows histograms of the clump masses for moderate 
resolution (left panel) and high resolution (right panel). At mod- 
erate resolution initially only a single clump forms, but seven or- 
bits later there are 5 bound clumps, all of several Ceres masses. 
The high-resolution run produces many more small clumps in 
the initial planetesimal formation burst. This is likely an eff'ect of 
the "hot start" initial condition where we turn on gravity during 
a concentration event as the high particle density allows smaller 
regions to undergo collapse. 



9. Summary and discussion 

In this paper we present (a) the first 512^ grid cell simulation 
of dust dynamics in turbulence driven by the magnetorotational 
instability and (b) a long time integration of the system per- 
formed at 256^ grid cells. Perhaps the most important finding 
is that large-scale pressure bumps and zonal flows in the gas ap- 
pear relatively independent of the resolution. The same is true 
for particle concentration in these pressure bumps. While sat- 
uration properties of MRI turbulen ce depend on the i i nplicit 
or explicit viscosity an d resist i yitv dLesur & Longarettii l2007l: 
iFromang & Papaloizoul l2007l: iDavis et al.L l2010l) . the emer- 
gence of large-scale zonal flows appears r elativel y independent 
of resolution (this work, Johansen et al., "2009a") and numeri- 
cal scheme (Fromang & Stone, 2009; Stone & Gardiner, 20T3). 
Particle concentration in pressure bumps can have profound in- 
fluence on particle coagulation by supplying an environment 
with h igh densities and low radial drift speeds (iBrauer et al.L 
'2008b), and on formation of planetesimals by gravitational con- 
traction and collapse of overdense regions (this work, J07). 

A direct comparison between the moderate-resolution and 
the high-resolution simulation is more difl&cult after self-gravity 
is turned on. The appearance of gravitationally bound clumps 
is inherently stochastic, as is the amplitude and phase of the first 
pressure bump to appear. The comparison is furthermore compli- 
cated by the extreme computational cost of the high-resolution 
simulation, which allowed us to evolve the system only for a few 
orbits after self-gravity is turned on. 

A significant improvement over J07 is that the moderate- 
resolution simulation could be run for much longer time. 
Therefore we were able to start self-gravity shortly after MRI 
turbulence had saturated, and to follow the system for more than 
ten orbits. In J07 self-gravity was turned on during a concen- 
tration event, precluding the concurrent evolution of self-gravity 
and concentration. This "hot start" may artificially increase the 
number of the forming planetesimals. The "cold start" initial 
conditions presented here lead to a more gradual formation of 
planetesimals over more than ten orbits. Still the most massive 
bound clump had grown to approximately 35 Mceres at the end 
of the simulation and was still gradually growing. 

The high-resolution simulation was given a "hot start", to 
focus computational resources on the most interesting part of 
the evolution. As expected these initial conditions allow a much 
higher number of smaller planetesimals to condense out of the 
flow. The most massive planetesimal at the end of the high- 
resolution simulation contained 7.5 Mceres of particles, but this 
"planetesimal" is accompanied by a number of bound clumps 
with masses from 0.5 to 4.5 Mceres- 
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The first clumps to condense out is on the order of a few 
Ceres masses in both the moderate-resolution simulation and 
the high-resolution simulation. The high-resolution simulation 
produced additionally a number of sub-Ceres-mass clumps. It 
thus appears that the higher-resolution simulation samples the 
initial clump function down to smaller masses. This observation 
strongly supports the use of simulations of even higher resolu- 
tion in order to study the broader spectrum of clump masses. 
Higher resolution will also allow studying simultaneously the 
concentration of smaller particles in sn ialler eddies and th eir role 
in the planetesimal formation process (ICuzzi et al l l2QTQh . 

We emphasize here the diff'erence between the Initial Mass 
Function of gravitationally bound clumps and of planetesimals. 
The first can be studied in the simulations that we present in 
this paper, while the actual masses and radii of planetesimal that 
form in the bound clumps will re quire the inclusion o f parti- 
cle shattering and part icle sticking. fZsom & DullemonJ (l2QQ8h 
and IZsom et al ] (1201011 used a representative particle approach 
to model interaction between superparticles in a 0-D particle-in- 
box approach, based on a compilation of laboratory results for 
the outcome of collisions depending on parti cle sizes, compo- 
sition and relative speed (iGiittler et al.l |201Q|) . We plan to im- 
plement this particle interaction scheme in the Pencil Code and 
perform simulations that include the concurrent action of hydro- 
dynamical clumping and particle interaction. This approach will 
ultimately be needed to determine whether each clump forms 
a single massive planetesimal or several planetesimals of lower 
mass. 

At both moderate and high resolution we observe the 
close approach and merging of gravitationally bound clumps. 
Concerns remain about whether these collisions are real, since 
our particle-mesh self-gravity algorithm prevents bound clumps 
from being smaller than a grid cell. Thus the collisional cross 
section is artificially large. Two observations nevertheless in- 
dicate that the collisions can be real: we observe planetesimal 
mergers at both moderate and high resolution and we see that the 
environment in which planetesimals merge is rich in unbound 
particles. Dynamical friction may thus play an important dissi- 
pative role in the dynamics and the merging. At high resolution 
we clearly see a trailing spiral arm exerting a negative torque on 
a planetesimal born in the vicinity of another planetesimal. 

If the transport of newly born planetesimals into each other's 
Hill spheres is physical (i.e. moderated by dynamical friction 
rather than artificial enlargement of planetesimals and numeri- 
cal viscosit y), then that can l ead to both mergers and production 
of binaries. iNesvorny et al.l (|201Q|) recently showed that gravi- 
tationally contracting clumps of particles can form wide sepa- 
ration binaries for a range of initial masses and clump rotations 
and that the properties of the binary orbits are in good agreement 
with observed Kuiper belt object binaries. 

In future simulations strongly bound clusters of particles 
should be turned into single gravitating sink particles, in order 
to prevent planetesimals from having artificially large sizes. In 
the present paper we decided to avoid using sink particles be- 
cause we wanted to evolve the system in its purest state with as 
few assumptions as possible. The disadvantage is that the par- 
ticle clumps become difficult to evolve numerically and hard to 
load balance. Using sink particles will thus also allow a longer 
time evolution of the system and use of proper friction times of 
large bodies. 

The measured a- value of MRI turbulence at 512^ is a ^ 
0.003. At a sound speed of Cs = 500 m/s, the expected colli- 
sion speed of marginally coupled m- sized boulders, based em- 
pirically on the measurements of J07, is ~ ^facs ~ 25 - 30 m/s. 



J07 showed that the actual collision speeds can be a factor of 
a few lower, because the particle layer damps MRI turbulence 
locally. In general boulder s are expecte d to shatter when they 
collide at 10 m/s or higher (lBenzll200Ql) . Much larger km- sized 
bodies are equally prone to fragmentation as random gravita- 
tional torques exerted by the turbulent gas exci te relative speeds 
higher t han the g ravitational escape speed (llda et al 1 |2008l 
Leinhar dt & Stewar t, 2009). A good environment for building 
planetesimals from bou lders may require a < 0.001, as in J07. 
I Johansen et al.l (l2009bh presented simulations with no MRI tur- 
bulence where turbule nce and particle cluniping is driven by the 
streaming instability (lYoudin & Goodman! l2005l) . They found 
typical collision speeds as low as a few meters per second. 

A second reason to prefer weak turbulence is the amount of 
mass available in the disc. If we apply our results to r = 5 AU, 
then our dimensionless gravity parameter corresponds to a gas 
column density of i7gas ~ ISOOgcm"^, mo re than ten times 
higher t han the Minimu m Mass Solar Nebula (I Weidenschillin gl 
Il977bl: iHayashil Il98ll) . Turbulence driven by streaming and 
Kelvin-Helmholtz instabilities can form planetesimals for col- 
umn densities compara ble to the Minimum Mass Solar Nebula 
(iJohansen et al.Ll2009bl) . 

The saturation of the magnetorotational instability is influ- 
enced by both the mean magnetic field and small scale dissipa- 
tion parameters, and the actual saturation level in protoplanetary 
discs is still unknown. Our results show that planetesimal forma- 
tion by clumping and self-gravity benefits overall from weaker 
MRI turbulence with a < 0.001. Future improvements in our 
understanding of protoplanetary disc turbulence will be needed 
to explore whether such a relatively low level of MRI turbulence 
is the case in the entire disc or only in smaller regions where the 
resistivity is high or the mean magnetic field is weak. 
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